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Abstract. The frozen Gaussian approximation, proposed in [Lu and Yang, 
|15|] . is an efficient computational tool for high frequency wave propagation. 
We continue in this paper the development of frozen Gaussian approximation. 
The frozen Gaussian approximation is extended to general linear strictly hy- 
perbolic systems. Eulerian methods based on frozen Gaussian approximation 
are developed to overcome the divergence problem of Lagrangian methods. 
The proposed Eulerian methods can also be used for the Herman-Kluk prop- 
agator in quantum mechanics. Numerical examples verify the performance of 
the proposed methods. 



1. Introduction 

This is the second of a series of papers on frozen Gaussian approximation for 
computing high frequency wave propagation. In the previous paper [15j . we pro- 
posed frozen Gaussian approximation for linear scalar wave equation with high 
frequency initial condition. It provides a valid approximate solution both in the 
presence of caustics and when the solution to wave propagation spreads. The frozen 
Gaussian approximation is based on asymptotic analysis in the phase space, and 
has better asymptotic accuracy than the Gaussian beam method [17] . The numeri- 
cal algorithm based on frozen Gaussian approximation was proposed in [15 within 
the Lagrangian framework. 

In the current paper, we provide an efficient methodology for computing high 
frequency wave propagation for general systems with smooth coefficients. On the 
one hand, we generalize the frozen Gaussian approximation to general linear strictly 
hyperbolic systems; on the other hand, we develop numerical methods based on 
the Eulerian formulation of frozen Gaussian approximation. The Eulerian methods 
solve the problem of divergence of particle trajectories in the Lagrangian method. 

Computation of wave propagation arises from many applications, for exam- 
ple seismology and electromagnetic radiation, where wave dynamics are governed 
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by hyperbolic equations. Direct numerical discretization of hyperbolic system 
is formidably expensive when waves are highly oscillatory. In conventional ap- 
proaches, the mesh size of discretization has to be comparable to wavelength or 
even smaller, while the domain of computation is determined by medium size. Dis- 
parity between the scales of wave length and medium size requires a huge number 
of grid points in each dimension. This makes computation extremely expensive. To 
bypass these difficulties of conventional approaches, numerical methods based on 
asymptotic analysis were developed, for example geometric optics and the Gauss- 
ian beam method. These methods are based on asymptotic analysis in the physical 
space, and solve an eikonal equation for the phase function S and a transport 
equation for the density p, 



where H(x,p) is the Hamiltonian function. The asymptotic expansion of geometric 
optics breaks down at caustics where the nonlinear eikonal equation (jl.ll) develops 
singularities. The Gaussian beam method replaces the real phase function S in 
geometric optics with a complex one, which makes asymptotic solution valid at 
caustics [20]. However, as discussed in [15], the construction of the Gaussian beam 
solution relies on Taylor expansion around beam center, therefore it loses accuracy 
when the beam spreads so that the width becomes large. 

Our previous work [15] made use of fixed-width Gaussian functions and carried 
out asymptotic analysis on phase plane to approximate the solution of high fre- 
quency wave propagation. It not only overcomes the shortcoming of the Gaussian 
beam method when beams spread, but also improves asymptotic accuracy. The 
numerical method given in [15j was of Lagrangian type, which may lose accuracy 
when particle trajectories are torn far way from each other after long time propa- 
gation. This divergence problem is of course a typical shortcoming of Lagrangian 
methods. One natural way to resolve it is to use Eulerian methods where numerical 
computation is done on fixed mesh grids. In the literature of high frequency wave 
computation, many Eulerian methods have been developed, for example wave front 
methods and moment-based methods reviewed in [1] , level set methods in geometric 
optics reviewed in [22] and the Eulerian Gaussian beam methods f5T4MlO|ITT j ll3 | ll4j . 
The underlying idea of these methods is to augment either geometric optics or the 
Gaussian beam method by numerical procedures based on partial differential equa- 
tions. In this paper, we propose Eulerian methods which augment frozen Gaussian 
approximation by numerical algorithms based on the Liouville equations, which is 
solved locally on the phase space. The proposed methods resolve the divergence 



(1.1) 
(1.2) 



d t S + H{x,V x S) =0, 
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problem in the Lagrangian method of frozen Gaussian approximation. As a byprod- 
uct, these Eulerian methods can be also applied to the Herman-Kluk propagator 
4 in quantum mechanics. 

The rest of the paper is organized as follows. In Section [2] we extend frozen 
Gaussian approximation (FGA) to general linear strictly hyperbolic systems. In 
Section [3l we discuss the application of FGA for high frequency wave propagation, 
emphasizing on the choices of parameters and discretization when the characteristic 
wave frequency is specified in initial conditions. The Eulerian formulation of frozen 
Gaussian approximation (EFGA) is introduced in Section[4] Two efficient numerical 
methods are proposed: Eulerian method and semi-Lagrangian method. Numerical 
results are shown in Section [5] We conclude with some remarks in Section [6] 

2. Frozen Gaussian approximation for general linear strictly 

hyperbolic systems 

We consider an M x M linear hyperbolic system in d dimensional space, 

d 

(2.1) d t u + J2Mv)d Xl u = Q, 

i=i 

where u = ( Ul ,...,u M ) T : R d -> R M and Ai : R d ->• M MxM , 1 < I < d are 
smooth M x M matrix valued functions in x. We assume that the system is strictly 
hyperbolic, i.e., for any p £ R d \{0} and any q G R d , the matrix Y^^iPl^lil) ^ as 
M distinguished eigenvalues, denoted as {H m (q,p)}f^ =1 . We denote by L m (q,p) 
and R m (q,p) the corresponding left and right eigenvectors, 

d 

(2-2) ^2 Pl Ll(q,p)Ai(q) = H m {q,p)L T m {q,P), 

1=1 

d 

(2.3) ^2piAi(q)R m (q,p) = H m (q,p)R m (q,p), 

i=i 

with the normalization 

lZ l {q,p)R n {q,p) = S mn , 

where 5 mn is the Kronecker delta function. As a result of the smoothness of A\ , the 
eigenvalues H m and eigenvectors L m , R m depend smoothly on (q,p). The method 
as presented requires only minor changes to be extended to hyperbolic system with 
eigenvalues of constant multiplicity; we will not go into details. 

2.1. Formulation. In frozen Gaussian approximation, to the leading order, the 
solution of the system (|2.1|) is approximated by the integral representation, 

i M r 

(2.4) u FGA (t,x) = d/2 / a m {t,q,p)e l ^ /£ v mfi {y,q,p)Aydpdq, 

^ ' m— 1 
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where v m ,o{y, q,p) = L m (q,p)uo(y) and i = \[—\ is the imaginary unit. Here we 
denote by Uq the initial condition of (|2.1j) . 
In (|2.4I) . the phase function $ m is given by 

(2.5) <f> m (t, x, y, q,p) = S m (t, q )P ) + l\ x - Q m \ 2 + P m ■ (x - Q m ) 

+ ^\y-q\ 2 p (y-q)- 

Here (Q m , P m ) are viewed as functions of t, q, and p. Given (q,p) as parameters, 
the evolution of (Q m (t,q,p), P m (t,q,p)) is given by the Hamiltonian flow with 
Hamiltonian function H m , 



(2.6) 



dt 

dP„ 

dt 



dp m H m (Q m , P m ), 

~ ®Qm {Q m , Pm ) , 



with initial conditions 

(2.7) Q m (0,9,p) -9, and P m (0,g,p)=p. 

The action function S m (t,q,p), also viewed as functions of (t,q,p), satisfies 

(2.8) ^ = P m ■ d Pm H m (Q m , P m ) - H m (Q m , P m ), 
with initial condition 

(2.9) S m (t,q,p)=0. 

The amplitude a m is given by a m (t, q,p) = a m (t, q,p)R m (Q m , P m ), where a m 
is determined by evolution equation (after dropping the subscript m for clarity and 
simplicity) , 

(2.10) ^ + aL T (d P H ■ d Q R-8 Q H ■ d P R) + a(d Zk L) T F.Z^ 1 

+ ad Zn Q 3 Z- r lL T (~d Qj A k + '-Pidg.dQ^R = 0, 

with initial condition 

(2.H) <r(P,q,p) =1 d ' 2 - 

In (|2.10p . we have used Einstein's summation convention and the short hand nota- 
tions 

(2.12) Fj = - ((Aj - d Pj H) + i{d Qj H - PidQ^)) R, 

(2.13) d z = d q -id p , Z = d z {Q + iP). 
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2.2. Asymptotic derivation. We justify the formulation of frozen Gaussian ap- 
proximation by asymptotics. We start with the following ansatz for the solution to 
(|2.1[) with the initial datum Uq, 

1 M f 

(2.14) u(t,x)= d/2 ^ (a m , (t,q,p) 

^ ' m — 1 

+ ea ma (t, g,p))e l * m/e -y m , (y, q,p) dy dpdq, 

where v m .o(y,q,p) — Lj n (q,p)uo(y), the phase function $ m is given in (|2.5p . and 
(Q m , P m ) follows the Hamiltonian flow (|2.6[) . 

We first state some lemmas that will be used later. The following lemma is 
essentially the same as that of Lemma 3.1 in |15j . and also the standard wave 
packet decomposition in disguise (see for example [2]). Hence the proof will be 
omitted. 

Lemma 2.1. For u e L 2 (R d ), it holds 

(2.15) u(x) = 1 / 2 d / 2 e^ ^y^u( y )dydpdq, 

where 

(2.16) $(0, x, y,q,p) = l -\x - q\ 2 + p ■ (x - q) + l -\y - q\ 2 - p ■ (y - q). 

Lemma 12.21 plays an important role in frozen Gaussian approximation. A similar 
observation in the case of the Schrbdinger equation was first noted by Kay [S] in 
asymptotic derivation of the Herman-Kluk propagator [4] in quantum mechanics. 
This observation was made precise in the work of Swart and Rousse [23] for the 
rigorous analysis of the Herman-Kluk propagator. Our previous work [15] extended 
it to linear wave equations. It is also true in the current case of general linear strictly 
hyperbolic systems. We omit the subscript m in the statement and proof of the 
lemma. 

Lemma 2.2. For any vector valued function b(y,q,p) and matrix valued function 
G{y 1 q,p) in Schwartz class viewed as functions of (y,q,p), we have 

(2-17) b(y,q,p) ■ (x - Q) ~ -ed^Z^ 1 ), 

and 

(2.18) (x-Q)-G(y,q,p)(x-Q) ~ e{d Zn Q ] )G ]k Z^+e 2 d Zr {d Zn {G 3k Z^u)ZJ r 1 ), 

where Einstein's summation convention has been used. 
Moreover, for multi-index a that \a\ > 3, 

(2.19) (x-Q) a -©(e 1 " 1 " 1 ). 
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Here we use the notation f ~ g to mean that 

(2.20) / fei* dydpdq = ge^Ay&pdq. 

Proof. Observe that at t = 0, 

d q S-(d q Q)P+p = 0, d p S - (d p Q)P = 0. 

Using (US]) and (HU), we have 

d t (d q S - (d q Q)P +p)=d q (d t S) - d q (d t Q) P - (d q Q)d t P 

= d q (P ■ d P H -H)- (d q (d P H))P + (d q Q)d Q H 
= {d q P)d P H - (d q Qd Q + d q Pd P )H + (8 q Q)d Q H 
= 0. 

Analogously we have d t (d p S — (d p Q)P) = 0. Therefore for all t > 0, 

(2.21) d q S-{d q Q)P + p = 0, d p S - {d p Q)P = 0. 

Then straightforward calculations yield 

d q <S> = (d q P - id q Q)(x - Q) - i{y - q), 
d p <Z> = {d p P - id p Q){x -Q)-(y- q), 

which implies that 

(2.22) id z <5> = Z(x - Q), 

where d z and Z are defined in (|2.13|) . The invertibility of Z follows the same 
argument in [151 Lemma 3.2], hence we omit the details here. 
Using ([2~2"2"|) . one has 

/ b- (x-Q)e^dydpdq^e [ bjZ^ 1 (-d Zk 3>) e** dy dpdq 

jR3d J R 3d \e / 

= -e[ (d^Zj^e^dydpdq, 

where the last equality is obtained from integration by parts. This proves (|2.17[) . 
Making use of (|2.17l) twice produces (|2.18l) 

(x - Q) ■ G(x -Q) = (x- ()K( - Q) k 

~-ed Zn ({x-Q) 3 G 3k Z^) 

= eid^Q^G^ - e(x - Q)jd z „(G jk Z^) 

~ eid^Q^Z^+e'd^d^kZ^Z^ 1 ). 

By induction it is easy to see that (|2.19p is true. 

□ 
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2.2.1. Initial value decomposition. We first check (|2.4|) gives the right initial value 
at time t = 0. Obviously it means to take, in (|2.14[) . 

a m ,o{0,q,p) = 2 d/2 R m (q,p) and a mA (0,q,p) = 0. 

We then have 

M 

«(0,x) = d/2 X! / 2 d / 2 J R m (g,p)iI(g,p)e^( ^'^^P)/ £ M o(y)dydpdq, 

^ ' m— 1 

where $(0, a;, y, q,p) is given in (|2.16[) . 

By the normalization of L m , R m , m = 1, . . . , A/, we have 

M 

7/7 = 1 

where Jm is the M x M identity matrix. Hence 

M(0,X) = (2^)3^/2 / 2 d/2 e lH °' X ' V ' q ' P)/£ My)<lydpdq = u (x). 

The last equality follows from Lemma |2. II Therefore the initial condition is repro- 
duced by El at t = 0. 



2.2.2. Evolution equation. We derive the evolution equation (|2.10j) for <r m in this 
subsection. Since the system under consideration is linear, we only need to consider 
one branch. For ease of notation, we suppress the subscript m in this section. 
Taking derivatives of $ with respect to t and x produces 

9t$ = d t S-P- d t Q + (x-Q)- (d t P - id t Q), 

and 

d Xl ® = i(xi - Qi) + Pi. 
Therefore the derivatives of ansatz ()2.14|) can be calculated as 

d t u = J (d t a +ed t a 1 + l -d t $(a Q + eai))e 4 * /e u (y, <?,£») dydpdq 

= J (~(d t S - P ■ d t Q)a + (d t a + i(d t S - P ■ d t Q) ai 

+ l -(x-Q)- [d t P - id t Q)a Q )y^/ e w(y, q,p) dydpdq + 0(e), 

and 



d Xl u = J (a a + eai)^-d Xl ^e i ' s ' /£ vo(y,q,p) dydpdq 



- £ Piao + {-\( x i ~ Qi)a + iPia 1 )^je l ' s ' /E vo(y,q,p) dydpdq + 0(e) 
Taylor expansion of Ai(x) around Ai(Q) gives 

A l (x)=A l (Q) + (x J -Q j )d Q] A l (Q) + ^(x ~Q ] )(x k ~Q k )d Q] d Qk A l (Q)+O(x-Q)\ 
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Substituting the above expressions into equation (|2.ip and matching orders in e 
yield the leading order equation, 

P d 

(2.23) / (d t S-P-d t Q + ^P l A l )a e l * /e v (y,q,p)dydpdq = Q. 
J 1=1 

Define the action function S to satisfy 

(2.24) d t S-P-d t Q = -H(Q,P), 
or equivalently 

(2.25) S(t, q,p)= f P-dtQ- H(Q, P) ds, 

Jo 

where P, Q and dtQ in the integrand are evaluated at (s, q,p). If we take 

(2.26) a (t,q,p) = a(t,q,p)R(Q,P), 
then by the definition of R(Q, P) in ((231) . 

d 

(d t S - P ■ d t Q + Y / PiAi)R{Q,P) = (d t S - P ■ d t Q + H(Q,P))r(Q,P) = 0. 

i=i 

To determine a, we investigate the next order equation, 



/('( 



dtS-P- d t Q + PiAijai + d t a 

+ ~(xj - Qj) (i{d t Pj - id t Q 3 )a Q - A 3 a + iPid Q] Aia^j 
(2.27) £ 

1 / i 

+ -(xj - Qj)(x k - Q k )[-d Qj A k a Q + -Pid Q] d Qk Aia Q 

xe l * /e w dydpdq = 0, 

where we interpret the terms quadratic in x — Q in the above expression by only 
keeping the 0(e) term arising from Lemma FZ7Z[ 

Solvability condition for a\ and Lemma 12.21 give the equation of ao, 
(2.28) 

L(Q, P) T (dta v Q (y, q,p) 

- d Zk (({idtPj + d t Qj)a ~ Aja + iPid Qj Aia )Zr k 1 v {y, q,p)j 

+ d Zn Qj(-d Qj A k aQ + ^Pid Qj d Qk Aia )Z~^v (y, q,p)j = 0. 

We next expand and simplify the above equation. For the first term, easy calcula- 
tions yield 

L(Q, P) T d t a = cV + aL T (dpH ■ d Q R - 8 Q H ■ d P R). 
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To simplify the second term in (|2.28p . we notice that by the definition (|2.3|) . 

d 

Y,PlMQMQ,P) = H{Q,P)R{Q,P). 

2=1 

Differentiating the above equation with respect to P and Q gives 

Ai(Q)R(Q,P) + PjAj(Q)d Pl R(Q,P) 

= d Pl H{Q, P)R{Q, P) + H{Q, P)d Pl R{Q, P), 

and 

P,d Ql A 3 (Q)R(Q 7 P) + P 3 A 3 {Q)d Ql R{Q, P) 

= 3 Ql H{Q, P)R{Q, P) + H(Q, P)d Ql R{Q, P). 
Taking inner product with L(Q, P) on the left produces 

(2.29) L T (Q 7 P){MQ) - d Pl H(Q, P))R(Q, P) = 0, 

(2.30) L T (Q, P)(P J 8 Ql A J (Q) - d Ql H(Q, P))R(Q, P) = 0. 

Recall the short hand notation 

Fj = (idtPj +d t Qj)R- A 3 R + iPid Qj AiR 
= -((A, - d Pj H) + i(d Qj H- Pid Qj Ai))R. 

Using (|2.29j) and (|2.30[) . it is clear that for any j = 1, . . . , d, 

L T {Q,P)F 3 =0. 

Hence, 

L T (Q,P)d Zk (aF 3 Z- k 1 v (y,q,p)) = aL T (Q, P)d Zk (Fj)ZJ k 1 v (y, q,p) 

+ L T (Q, P)F,d Zk (aZ^voiy, q,p)) 
= -*(d Zk L) T (Q, P)F j Zr k 1 v (y, q,p). 

Therefore, (|2.28p can be rewritten as 

d t <7 + aL T (d P H ■ OqR - d Q H ■ d P R) + a(d Zk L) T F ^Zj^ 

+ ^ Zn QjZ^L T {-d Qj A k + '-P^Oq^R = 0, 

which is just the evolution equation (|2. 10[) . 

2.3. Examples. We apply the general results to some specific systems. For a given 
system, once the eigenvalues and eigenfunctions are determined, it is straightfor- 
ward to obtain the initial value decomposition and evolution equation for a. We 
illustrate this by two examples. 
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2.3.1. Scalar wave equation in one dimension. Consider the ID scalar wave equa- 
tion 

(2.31) d 2 u- c 2 {x)d 2 x u = 0, 

where c(x) > is the (local) wave speed. Define r = dtu and s = d x u, and 
transform (|2.31[) into the system 



d t r - c{x) 2 d x s = 0, 
dts — d x r = 0. 



(2.32) 

It can be rewritten as 
where A is given by 

The eigenvalues of A are given by 

H±(q,p)=±c(q)\p\. 

Hence the 2x2 system (I2.32|) is strictly hyperbolic. The corresponding right and 
left eigenvectors are 

By (|2.10|) , the evolution equations are given by 

dta± = ±^^ C 'W±) ± ^-Z±d*Q± (2^c'(g ± ) - i\P ± \c"(Q ± )) . 

This agrees with the evolution equations given in [15] in one dimensional case, 
where the amplitude was denoted as a instead of a. 
In (|2.4I) . the initial value decomposition is taken as 

Vo,±(y,q,p) = jr( , i , . d t u{0,y) T -d y u(0,y)). 

2\\p\c{q) p J 

If the initial value to the wave equation takes the WKB form, i.e., 

u (x) = A a (x)ei s °( x \ 
d t u (x) = ±B (x)ei s °(*\ 



then 



«D,±(»,?,P) = ^(^ e!S0b)/£ ^( !A *) S »fe)+^(!')) e,S0(j)/e 

Remark. The choices of R and L are not unique. The above choice is made in 
order to match the results in [15]. If different normalization is chosen for R, the 
results of initial value decomposition and evolution equations can be different. 
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2.3.2. Acoustic wave equation in two dimension. We next consider the acoustic 
wave equation in two dimension, 



(2.33) 



d t V + VU = 0, 
d t n + c 2 (x)\7 ■ V = 0, 

where V is velocity and II is pressure. Define u = (Vi, V<x, n) T , and we can rewrite 
(|2.33p as a 3 x 3 linear hyperbolic system, 

d t u + A 1 (x)d Xl u + A 2 (x)d X2 u = 0, 

where 





( o 









fo 















3 


A 2 = 








3 




\c{xf 











c(x) 2 





Then the eigenfunctions in (|2.2f) - (|2.3|) are given by 

H 1 (q,p) = 0, H ± (q,p)=±c(q)\p\, 

where we have used the subscripts ± instead of number subscripts. This implies 
the system is strictly hyperbolic, and the corresponding eigenvectors are 

±Pi 

±Pl 

L / ±pi/\p\ 2 

I ±Pl/\p\ 2 

' \l/(c(q)\p\), 



Li(q,p) = r-|2 



/ Pi 






-PI 




V 






1 


( 


Pi 








-Pi 


w 


{ 








L±(q,p) 



In (|2.4I) , the initial value decomposition is taken as 

«o,i(v,g,p) = T^{piV 1 °(y)-piV 2 °(y)), 



\Pl 



and 



«0,±(V,9,P) = ^(±PiV 1 "(y)±PiV 2 (y) + ^ II (</)), 

where Uq — (V®, V 2 , II ) T is the initial condition. 

The evolution equation (|2. 10[) of a can be simplified as, after straightforward but 
lengthy calculations, 

da i 



dt 



0. 



and 



d<7± 
~dT 



<r± ( P± 



o-± 



±^tv[Z ± 'd z Q ± [2 



ic 

\pT\ 
p± 



\p±, 
p±®p± 
\p±\ 2 



d Q± c 



l)-i\P ± \d' Q± c 
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We note that the solution associated with the first branch (Hi — 0) does not involve 
in time, since the Hamiltonian flow (|2.6|) is Qi(t) = q, Pi{t) = p and a\ stays 
constant. 

3. High frequency wave propagation 

The frozen Gaussian approximation (FGA) formulated in Section[3]approximates 
the propagation operator of hyperbolic system. The approximation is useful espe- 
cially in the case of high frequency wave propagation, where the small parameter e 
should be chosen according to the initial condition. 

We rewrite the frozen Gaussian approximation (|2.4I) as 

1 M f 

u(t,x) = (27r , 3d/2 / a o^m{t,q,py' 5 ' m/e Ll l (q 1 p)u (y)dydpAq 

* ' m— 1 

l M r 

= {2lT£ fd/2 E J «o,m(*,Q,p)e 8 



mi i M 

'" I / , ^ Sm +i\a : -Q m \' 2 /2+P rn -( x -Q m ))/s 



x Lj l (q,p)ilJ (q,p)dpdq, 

with ipo(q,P) given by 

(3.2) *l> Q (q,p) = I e-\*-rf/W-^y-«V s u (y)dy. 



This implies that the initial condition is first transformed to be on the (q,p) 
phase plane by taking inner product with Gaussian functions, then one evolves the 
centers (Q m , P m ) and weights a m of these Gaussian functions by p.6p and (|2.10|l . 
At final time t, the solution is approximated by superpositions of these Gaussian 
functions. We address two issues appearing in numerical algorithms of FGA: one 
is to estimate the mesh size of (q,p) in the discretization of Q3.1J : the other is to 
discuss the choice of small parameter e when the initial condition takes the form 

(3.3) u(0,x) = An^Je 15 "^)^, 

where Aq and So are smooth and compact support functions. The parameter r\ 
characterizes the frequency of the initial wave. Small rj indicates high frequency 
waves. 

3.1. Mesh size of (q,p). Taking derivatives of (|3.2|) with respect to p and q 
produces 

(3.4) 9 p ^ (q,p) = — <p(q,p); 

1 1 

(3.5) d q ip (q,p) = -ip(q,p) + -p®i/> (q,p), 



where 



f(q,p) = J(y~q)® e-\y-^/^-^-^ e u {y) dy. 
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Denote the integrand in (|3.ip as, where we drop the subscript m without loss of 
generality, 

(3.6) A(t,x,q,p) = a (t, q ,p)e^ S+ ^-^ 2 / 2+P ^- Q ^/ £ L T (q,p)xP Q (q,p). 
Taking derivatives of (|3.6[) with respect to p and q yields 
P A = a p ai T ^ e l(S+l|a! - Q|2/2+jP - (a! - Q))/E 

+ (~ i d p s + ld pQ -(Q-x)+ d p P ■ (x - Q) - d p Q ■P)L-ip 

+ d p L ■ Vo - % -<P ■ Lj ® ae <S+A X -Q?/-2+P-( X -Q))/e^ 

and 

d q A = d q aL T ^ s +^-^ 2 / 2+p < x -^/" 

+ {^( d i S + ld iQ -(Q- x ) + d g P -{x-Q)-d q Q-P)L- O 

+ d q L ■ -0o + \V ■ L + -pL ■ O ) ® ae <S + r\ X -Q\V2+P-( X -Q))/e^ 

By (|2.2ip in the proof of Lemma 12.21 we can simplify the above expressions of 
derivatives as 

d p A = dpaL T iP e< s +^-W 2 / 2 + p < x -Wc 

-(d p P — idpQ) ■ (x — Q)L ■ O 

+ dpL ■ll, Q - % -<p-L S j® ae <S + A X -Q?/2 + P-( x -Q))/s^ 

and 

d q A = d q aL T ^e< s +^-^ 2 ^+ p < x -^/ £ 
-(d q P-id q Q) ■ (x-Q)L-^ 
+ 8 q L • O + ~v ■ L\® ae <s+A*-Q\ 2 /2+P-(*-Q))/e^ 
Keeping only the highest order terms gives 

d p A = (^-(dpP - idpQ) ■ (x - Q)L ■ O - -<p ■ Lj ® a 

x e i(S+i\ X -Q\ 2 /2+P-( x -Q))/s + 0(1), 

and 

d q A = (± ({d q P - id q Q) ■ (x - Q)L • O + ■ Lj ® a 

x e i(S+i\ X -Q\ 2 /2+P-( x -Q))/e + Qn\ 

Notice that (x — Q) and (y — q) is 0(e 1 ^ 2 ) due to the Gaussian factor, therefore 
both derivatives are 0(e~ 1 ^ 2 ), while the function A is 0(1). As a result, in order 



e 



14 



JIANFENG LU AND XU YANG 



to get an accurate discretization of the integral (|3.1jl , one has to take the mesh size 
in q and p to be at least ©(e 1 / 2 ). 

3.2. Choice of parameter e. While the original hyperbolic system lives in phys- 
ical domain, FGA works on phase plane, hence the dimensionality is doubled. The 
cost of numerical algorithm based on FGA can be estimated by the number of mesh 
points used on phase plane. This means we need to find the region where xp (q,p) 
makes a significant contribution. In this subsection we investigate the effect of e 
on the size of region under the consideration of the high frequency initial condition 
(pLl. 

Substitute the initial condition (|3.3|) into (|3.2[) . we have 



^o(q,P) = J e-^-^ 1 7(2e)-«P-(»-9)/«+»So(y)/u A {y)dy. 

We will choose e comparable to r\ and discuss the effects of increasing or decreasing 
e. Let r = e/rj, then 

ip (q,p) = J e-\ y - q \ 2/{2s] e l( - p - {y - q)+rSo{v))/e A {y)&y. 

Taylor expansion gives 

So(y) = S (q) + V5 (q) ■ (y - q) + \^ 2 S (q + 9(y - q)) : (y - qf , 

where 8 € [0,1] depends on y. Define 

Mv, q) = s Q (y) - s ( g ) - vs (q) ■ (y - q), 

then we have 

ip (q,p) = e wSa ^' E J e-l«-JI'7(^) e '(-P+ I ' vs ''('?))-(«- < J)/ £ e lflo( ^ <3)/£ A ( ? /)d 2 / 



el H, ( q +Vty, q )/e AQ ( q + dy 



Define 



f q (y) = e -i«i 3 /2 e ^o( g +^, q )/ £j4o( g + ^ y) _ 

By the definition of R a , it is clear that the derivative of / with respect to y is 
bounded independent of e. Therefore, by 

(3.7) ^ (q,p) = V~ee trS ^f q (e-^ 2 (p-rWS (q))), 

standard integration by parts argument yields 

(3-8) \MQ,P)\ < ( 



\p-rVS (q[ 



IJV' 



_(JV+l)/2 



for any positive integer N. In (|3.7|) . f q means the Fourier transform of f q , and 
(|3.8p is actually the decay rate of the Fourier transform. 
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The equation (|3.8j) implies xp (q,p) makes significant contributions only when 
\p — rVSo{q)\ is 0(e 1 ^ 2 ). Therefore one only needs to consider the region where 
p is localized around rVSo(<z). If one takes the mesh size of p as 0(e 1 ^ 2 ), then 
the number of mesh points in p given q is a constant. The total number of (q, p) 
points to be considered is 0{e~ d / 2 ) ~ 0(rj~ d / 2 ). Therefore while smaller e gives 
better asymptotic accuracy, it requires more computation cost. This is a trade-off 
between cost and performance. 

It is seen that FGA is suitable for computing high frequency wave propagation 
when e is taken comparable to 77, which is the reciprocal of the characteristic fre- 
quency of initial wave field. We remark that, in a follow-up work j!6j . we establish 
rigorous analysis on the accuracy of FGA for high frequency wave propagation for 
general linear strictly hyperbolic system. 

4. Eulerian Frozen Gaussian approximation 

In this section we introduce Eulerian frozen Gaussian approximation (EFGA) 
for computation of linear strictly hyperbolic system. We first describe Eulerian 
formulation, followed by numerical algorithms based on the Eulerian formulation. 
These Eulerian methods can also be applied for computation of the Herman-Kluk 
propagator [4] in quantum mechanics, which is discussed in the last subsection. 

4.1. Eulerian formulation. The formulation of EFGA is given by 

M 

(4.1) u EFGA (t,x) = —^ m <r m (t,Q,P)R m (Q,P)e^dPdQ, 
where the phase function B m is 

(4.2) e m (t,x,Q,P) =S m (t,Q,P) + P- {x-Q) + ~\x-Q\ 2 . 
Define the Liouville operator 

£-m = dt + dpH m ■ 8q — dqH m ■ dp. 
The evolution of S m (t, Q, P) satisfies 

(4.3) C m S m = P ■ dpH m — H rn , 

with initial condition S^O, Q, P) = 0. 

To get the evolution of a m , we define the auxiliary functions 

<l> m (t,Q,P) = (0m,l,-" ,<f>m,d), 

given by 

(4.4) £ m m = 0, 
with initial condition 



(4.5) 



<j> m {0,Q,P) = P + iQ. 
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Once (f) m is determined, the evolution of a m (t,Q, P) is given by (where we have 
omitted the subscript to for simplicity of notation) 

Ca = - crL T (dpH ■ d Q R - d Q H ■ d P R) 

(4.6) - ^{dp<t>k ■ d Q L - d Q <j) k ■ d P L) T FjZ7* 

- ad Pi cfi n Z^L T (-d Qj A k + ^PidQ j dQ k Ai)R, 

with initial condition 

(4.7) a m (0,Q,P)=2 d / 2 J v mfi (y,Q,P)exp(^(-P-(y-Q)+ t -\y-Q\ 2 )^dy, 

where u m ,o(2/> Q> P) = Lj n (Q , P)u a (y) . In (|4.6[) . we have used the shorthand 
notations, 

(4.8) Z= (9p0) T - l (a Q 0) T , 

(4.9) Fj = - ((A* - d Pj H) + i{d Qj H - PidQjAt)) R. 

4.2. Derivation. Under the change of variable, 

R d x R d — > R d x R d 

(q,p) — > {Q,P) 

where r) m is the Hamiltonian flow given by (I2.6[l - (|2.7p . the FGA formulation (12. 4p 
can be rewritten as 

M 

u EFGA (t,x)= —— / a m (t,Q,P)e^^v m ^y,q m ,p m )dydPdQ, 

^ ' m—l 

where (q m ,p m ) = ?7~ 1 (<5, P) and we have used the fact that the Jacobian of r] m 
equals to 1 due to symplecticity of the Hamiltonian flow. 
The phase function $ m is given by 

(4.11) $ ro (t, x, y, Q, P) = S m (t, Q, P) + l -\x - Q\ 2 + P ■ (x - Q) 

+ \\y-<i m \ 2 -p m ■ (v-Qm)- 

Remark. In the Lagrangian formulation of FGA, solution of each branch starts at 
the same (q,p), so that (q,p) is independent of m while (Q,P) given by (|2.6p 
depends on to; in the Eulerian formulation of FGA, solution of each branch ends 
at the same (Q, P), therefore (Q, P) is independent of to while (q, p) given by Ty" 1 
depends on to. 

What remains is to derive evolution equations of a m and S m in terms of Eulerian 
coordinates (Q, P). The easy observation is to change time derivative in Lagrangian 
coordinate to the Liouville operator in Eulerian coordinate by chain rule and 

-77 — ► £m = dt + dpH m ■ 8q — dQH m ■ dp, 
which, together with (|2.8[) . implies (|4.3|) . 
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The difficult part is, in the evolution equation (|2.10[) for a m , there are terms 
containing Lagrangian derivatives with respect to q and p. To replace them in 
Eulerian coordinate, one needs the following theorem, where we omit the subscript 
m for simplicity. 

Theorem 4.1. Assume 0(t, Q, P) = (<pi, • • • , cj)^) is the solution to 

(4.12) d t <l> + dpH d Q( p-d Q H ■d P <f> = 0, 

with initial condition 

(4.13) 0(0, Q,P) =P + iQ. 
Denote 

X = d z Q, Y = d z P, 
where d z Q and d z P are given in Lagrangian coordinate, then 

(4.14) X=(d P <f>) T , Y = -(d Q cj>) T , 



in 



Eulerian coordinate. 



Proof. Differentiating (|2.6[) with respect to z produces 
(4.15) 



dX v d 2 H ^8 2 H 
= X——— + Y 



dt dQdP dP 



2 1 



dY _ d 2 H d 2 H 
~dT ~~ X ~dCf~ dPdQ 1 
which are the equations for X and Y in Lagrangian coordinate. 
Therefore X and Y satisfy, in Eulerian coordinate, 



(4.16) 



2 ! 



CX - X dQdP +Y W 



dQ 2 ' dPdQ' 
Differentiating (|4.12|) with respect to P and Q yields 



(4.17) 



C ^=-dQdP dQ<j> + W ^ 



Comparing (|4.16j) with ()4.17j) . one can see that X and Y satisfy the same the 
Liouville equations as (dp<f>) T and — (9q</>) T . Moreover, (|4.13[) implies that X and 
Y have the same initial conditions as (dp<p) T and — (9q</>) T . Since the Liouville 
equation is linear, we have in Eulerian coordinate, 



x = {dp ( t>) T 1 Y = -(d Q 4>) T . 



a 
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Remark. We observe that (|4.15j) is equivalent to the dynamic ray tracing equations 
in Gaussian beam method [17j . Therefore Theorem l4.1l can be applied in computing 
Hessian functions in Eulerian Gaussian beam methods, which is essentially the 
approach used in [5HH] . 

Hence (l2~T0)) and (|4T4| imply (|4~6)l . 
Rewrite $ m in (|4.11|) as 

$m = e m + -\y -q m \ 2 -p m ■ (y-q m ), 

where 8 m is given by (I4.2[) . Since q and p are parameters in (|2.10l) . the initial 
condition (|4.7[) is obtained by combining the term 

J v m , (y, q,p)e-^<v-9nJ-h\v-i m \ a dy 
with the initial condition (|2.11[) . 

4.3. Algorithm. We introduce two numerical algorithms based on the Eulerian 
formulation: Eulerian method and semi-Lagrangian method. Let us first describe 
the meshes needed in these algorithms. 

4.3.1. Numerical meshes. 

(1) Discrete meshes of Q and P for solving the Liouville equations. 

Denote SQ = (SQi, ■ ■ • , SQd) and SP = (6 Pi, ■ ■ ■ , SPd) as the mesh size. 
Suppose Q° = (Qi, • • • , Q d ) is the starting point, then the mesh grids Q k , 
k = (ki, ■ ■ ■ , kd), are defined as 

Q k = {Qi + (ki - l)SQi, ■■■ ,Q° d + {kd~ l)5Q d ), 

where kj = !,■■■ ,N q for each j £ {1 , • • • , d} . 

The mesh grids P e , £ = (£ u • ■ • Jd), are defined as 

P* = (A° + (h - l)SPi, ...,pg + (i d - l)6P d ) , 

where £j = 1, • • • ,N P for each j £ {1, • • • , d} and P° = (Pf, ■ • • , P°) is the 
starting point. 

(2) Discrete mesh of y for evaluating the initial condition in (|4.7[) . Sy = 
(Syi, • • • , 8yd) is the mesh size. Denote y° = • • • , y d ) as the starting 
point. The mesh grids y m are, m — (mi, ■ ■ ■ , rrid), 

y m =(yi + (m 1 -l)Sy 1 ,--- ,y°d + (m d -l)5yd), 

where mj = 1, ■ ■ ■ , N y for each j £ {1, ■ ■ ■ , d}. 

(3) Discrete mesh of x for reconstructing the final solution. Sx = (6xi , Sxd) 
is the mesh size. Denote x° = (x®, ■ ■ ■ , x d ) as the starting point. The mesh 
grids x n are, n — (m, ■ ■ ■ , rid), 

x n = (a;? + (ni - l)Sxi,--- ,x° d + (n d - l)5x d ), 
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where rij = 1, • • • , N x for each j € {1, • • • , d}. 

4.3.2. Eulerian method. With the prepared meshes, the Eulerian frozen Gaussian 
approximation algorithm is given as follows. 

Step 1. Compute L m (Q k , P £ ) and R m (Q k , P e ) by solving the eigenvalue problems 
(l2~2l- (l2~3f . 

Step 2. Compute the initial condition (|4.7|) at (Q k ,P £ ), 

(4.18) a m (0, Q k , P e ) = 2 d ' 2 ]T e i(^-(y m -Q k )+i\v m -Q k \ 2 ) 

m 

x v m , (y m , Q k ,P g )re(\y m - Q fc |)<5yi • • • Sy d , 

where rg is a cutoff function such that rg = 1 in the ball of radius > 

centered at origin and rg = outside the ball. 
Step 3. Solve (|4.3[) . (|4.4[> and (14.61) by finite difference/ volume/element methods, 

for example, standard upwind scheme with van Leer flux limiter f |12|b 

Denote the numerical solutions as S%> and o-^f. 
Step 4. Reconstruct the solution by (|4.1|) . 

(4.19) i=it^ 27r£ ) 7 

x fi m (Q fe , P'Jrfldx" - Q fc |)<5g 1 • • • SQ^Pi • • • SP d . 

Note that a naive implementation of the above method will result in numeri- 
cal methods on the phase plane, and hence doubles the dimensionality. A more 
efficient way is to implement this method locally on the phase plane, when initial 
conditions have localization properties. This local solver strategy is important to 
make Eulerian methods efficient, and we detail the algorithm below. 

If one considers WKB initial condition for linear strictly hyperbolic system (|2.1I) , 



(4.20) u(0,x) = a (x)e 



) 

then the initial condition (|4.7[) is localized in momentum space around the sub- 
manifold P = VqS'o(Q) as discussed in Section I3~2"1 A one-dimensional example is 
given in Figure [TJ This localization property allows efficient local implementation 
of the Eulerian method. One possible and simple strategy is based on indicator 
functions. The idea is similar to the moving mesh algorithm [3]. 
Define indicator functions n m , m = 1, • • • , M, which satisfy 

(4.21) £ m K m = 0, 

with initial condition 



(4.22) Km (0,Q,P) 



1 ifa m (0,Q,P) ^0, 
otherwise. 
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0.5 1 1.5 



Q 

Figure 1. An illustration of the localization of <r(0, Q, P) in one 
dimension for v (y) — & exp , £ — 1/128; the black solid 

curve is P — cos(6Q)/2. 

Then in Step 3 of the algorithm, when solving (|4.3I) , (14.41) and (I4.6[) one only needs 
to update function values on those (Q k , P ) where n m is nonzero. 

Remark. 1. In setting up the meshes, we assume that initial condition either has 
compact support or decays sufficiently fast to zero as x —¥ 00 so that we only need 
finite number of mesh points in space. 

2. The role of the truncation function r e is to save computational cost, since 
although Gaussian function is not localized, it decays quickly away from the center. 
In practice we take 8 — 0(y/F), the same order as the width of each Gaussian, when 
evaluate (|4.7[) and (|4.1I) numerically. 

3. There are two types of errors present in the method. The first type comes 
from the asymptotic approximation to strictly hyperbolic system. This error can 
not be reduced unless one includes higher order asymptotic corrections. The other 
type is the numerical error which comes from two sources: one is from solving the 
Liouville equations numerically; the other is from the discrete approximation of 
integrals (|4.7|) and (|4.1j) . It can be reduced by either taking small mesh size and 
time step or using higher order numerical methods. 

4. Step 2 and 4 can be expedited by making use of discrete fast Gaussian 
transform, as in [T8lll9| . 

4.3.3. Semi-Lagrangian method. Alternatively, one can also use semi-Lagrangian 
method. This is a type of Lagrangian method based on the Liouville equations (|4.3|) . 
(|4.4p and (|4.6|) , which can be viewed as a local implementation of Eulerian method 
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on an adaptive mesh. Different from the meshes of Eulcrian method in Section 
14.3-H the meshes of (Q, P) in semi-Lagrangian method is determined adaptively 
from initial conditions, while the meshes of y and x are still the same. 

The underlying idea is to first lay down uniform mesh grids of (Q,P), evolve 
the grids to time t according to Hamiltonian flow (denoted by (C^P 4 )); then set 
up uniform mesh grids of (Q,P) at time t based on (Q*,P*) and use method of 
characteristics to compute the solutions to the Liouville equations. The difference 
from Eulerian method is that it solves the Liouville equations by numerical inte- 
grators for ODE instead of numerical schemes for PDE. The detailed algorithm is 
given as follows, where we only focus on computation of one eigenvalue branch and 
omit the subscript m for simplicity and clarity. 

Step 1. Choose initial uniform mesh grids (Q s '° , P r '°) where er(0, Q 
nonzero. Solve time-forward Hamiltonian flow 

(^ = d P H(Q,P), 
^ = -d Q H(Q ) P), 



s,0 pr.0) ig 



with initial conditions 



Q(0,Q s ' ,P r ' ) 



Q s 



and P(0, Q s 



jr.O 



Denote (Q^P 7 ^) = (Q(t, Q s '°, P r '°), P(t, Q s '°, P r ' )) . 
Step 2. Choose uniform mesh grids (Q k ,P £ ) so that all the points (Q s '*,P r,t ) lie 
in mesh cells. Solve time-backward Hamiltonian flow 

( ^- = -d P H(P,Q), 
at 

¥L=d Q H(P,Q), 

with initial conditions 

Q(0,Q k ,P e ) = Q k , and P(0, Q k , P e ) = P £ . 
In the meantime, solve time-backward equation 

-P-dpH + H, 

with initial condition S^O, Q fc , P £ ) = 0, then 

S(t,Q k ,P £ ) = -S(t 1 Q k ,P e ). 
Denote (q k ,p £ ) = (Q(t,Q k ,P t ),P(t,Q k ,P i )), then 
<j>(t,Q k 1 P e )=p e + tq k . 



~dt 
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Step 3. Solve time- forward equation 

4? + aL T (d P H ■ d Q R - d Q H ■ d P R) + a(d Zk L) T F.Z. 1 
at 

+ od Zn Q,Z-^l7(-d Q] A k + "-Pidg.dQ^R = 0, 

where 

F 3 = - ((A 3 - d Pj H) + i(d Qj H - Pidq^j) R, 
with initial condition 

*(0,q k ,p*) = 2 d / 2 J v (y,q k ,p e ) 

x exp f 1 ( - p e ■ (y - q k ) + l -\y - q k \ 2 )^j dy. 

Then a(t,Q k ,P £ ) = a(t,q k ,p £ ). 
Step 4. Reconstruct the solution by 

x R(Q k , P*)re(\x n - Q k \)SQ 1 ■ ■ ■ SQdSPx ■ ■ ■ SP d , 
where R the right eigenfunction in (|2.3[) . 

Remark. In Step 3 one needs to compute Z = d z (Q + iP). Since (q k ,p £ ) is not 
a uniform mesh grid, it can lose accuracy by using divided difference to compute 
d z Q and d z P. To resolve this problem, one can solve f)4. 15[) to get d z Q and d z P 
instead. 

4.4. Eulerian method for the Herman-Kluk propagator of the Schrodinger 
equation. The Eulerian method introduced in Section 14.11 can be also applied to 
computation of the Herman-Kluk propagator in quantum mechanics. 
The rescaled linear Schrodinger equation is given by 

d^ £ e 2 j 
(4.23) ie— = - — A^ e + U{x)^ e , xeR d , 

where ^ e (t,x) is the wave function, U(x) is the potential and e is the re-scaled 
Plank constant that describes the ratio between quantum time/space scale and the 
macroscopic time/space scale. This scaling corresponds to the semiclassical regime. 

We briefly describe the formulation of the Herman-Kluk propagator [3] below. 
One can see that it is similar to the frozen Gaussian approximation. In fact, the 
frozen Gaussian approximation introduced in [15] is motivated by the ideas of the 
Herman-Kluk propagator. We remark that semiclassical approximation underlying 
the Herman-Kluk propagator has been recently rigorously analyzed by Swart and 
Rousse [23] and Robert [2"T] . 
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Using the Herman-Kluk propagator, the solution to the Schrodinger equation is 
approximated by 

(4.24) * £ HK (*.*) = (27re L/a / a(t,q,p)e^y^^%(y)dydpdq, 
where is the initial condition of (|4.23l) . Here the phase function $ is given by 

(4.25) $(t,x,y,q,p) = S{t,q,p) + l -\x - Q\ 2 + P ■ (x - Q) 



^\y-q\ 2 -p(y-q)- 



The evolution of (Q, P) satisfies 
(4.26) 




with initial conditions 

(4.27) Q(0,q,p) = q, and P(0,q,p)=p. 

The action function S(t,q,p) satisfies 

with initial condition 5(0, q,p) = 0. The evolution of a(t, q,p) satisfies 

(4.29) = ~<ztr (z-\d x P - id z Qd 2 Q U) 

with initial condition a(0, q,p) = 2 d l 2 . 

It is easy to see that (|4.26[) and (|4.28p are the same as (|2.6[) and (|2 . 8[) if one 
takes H m (Q,P) = \P\ 2 /2 + U(Q). The difference lies in the amplitude evolution 
equation (|4.29p . But it does not raise any difficulty for Eulerian formulation. One 
can still write down Eulerian formulation based on Theorem 14. 11 

(4.30) *| HK (i, x) = {2 J )3d/2 J o(t, Q, P)e i0 / £ J- 1 dPdQ, 
where the phase function G is 

(4.31) 0(t, x, Q, P) = S(t, Q,P) + P ■ (x - Q) + % -\x - Q\ 2 . 
Define the Liouville operator 

C = d t + P-d Q -d Q U- Op. 
Then the evolution of 5(t, Q, P) satisfies 

\P\ 2 

(4.32) CS = J-L - U{Q), 
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with initial condition S(0, Q, P) — 0. We introduce the auxiliary function <p(t, Q, P), 
which satisfies 

(4.33) Ccj> = 0, 
with initial condition 

(4.34) 0(0, Q,P) =P + iQ. 
With <p determined, the evolution of a(t, Q, P) satisfies 

(4.35) £a = -^atr [Z' 1 ((d Q <f>) T + i{d P ct>) T d 2 Q U)^j , 
where 

(4.36) Z = (dpc^f - t(d Q( j)) T . 
The initial condition of (|4.35[) is prepared as 

(4.37) o(0, Q, P) = 2 d ' 2 J %(y) exp f~( - P ■ (y - Q) + % -\y - Q| 2 )) dy. 



Therefore the numerical methods discussed in section 14731 can also be applied to 
the Herman-Kluk propagator. 

5. Numerical examples 

In this section, we present four numerical examples to show the performance of 
Eulerian frozen Gaussian approximation (EFGA) and also Eulerian methods for 
the Herman-Kluk propagator. Two of the examples correspond to EFGA of wave 
propagation discussed in Section 12.31 and the other corresponds to the Schrodinger 
equation discussed in Section 14.41 We consider WKB initial conditions, and use 
the Eulerian method with local indicator to compute wave propagation, and use 
the semi-Lagrangian method for the Herman-Kluk propagator of the Schrodinger 
equation. 

5.1. Wave propagation. 

Example 5.1 (One-dimensional scalar wave equation). 

d 2 t u - c{xfd 2 x u = 0. 
The wave speed is c(x) — x 2 . The initial conditions are 

uo = exp (-100(a; - 0.5) 2 ) exp (— ) , 

d t u = exp(-100(x - 0.5) 2 ) exp — 

£ V £ 
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Table 1. Example ED the £°° and i 2 errors for EFGA. 



£ 


1/2 6 


1/2 7 


1/2 8 


|| U _ U BFGA||^ 


, 1.46 x 10" 1 


6.39 x 1(T 2 


3.52 x 10~ 2 


\\u-u™ GA \\ e 


4.81 x 10~ 2 


2.39 x 1(T 2 


9.11 x icr 3 



The final time is T = 0.8. We plot the real part of the wave field obtained by 
EFGA compared with the true solution in Figure El for e = 1/64, 1/128, 1/256. As 
one can see, the span of the solution reaches 2 roughly at T — 0.8, although it starts 
with only 0.5 approximately. Apparently the wave spreads quickly in this example. 
Table Q] shows the l°° and I 2 errors of the EFGA solution. The convergence orders 
in e of £°° and I 2 norms are 1.02 and 1.20 separately for EFGA, which confirms the 
asymptotic accuracy. The true solution is computed by the finite difference method 
using the mesh size of N x — 49152 and the time step of N t = 524288 for domain 
[0, 6]. We take N t = 1024, Sq = Sp = Sy = 1/128 and 5x = 1/2048 in EFGA. 

Example 5.2 (Two-dimensional acoustic wave equations). 

jd t V + VII = 0; 
[d t n + c 2 (x)V ■ V = 0, 

where V — (Vi, V2), x = (x\,X2) and c(x) = 1. The initial conditions are 

n = \J 1 + sin 2 (4x 2 )/16exp(-100(^ + x\)) exp Q( - x t + cos(4x 2 )/16)) , 
Vj.,0 = -exp(-100(x 2 + x 2 )) exp Q( - x x + cos(4ar 2 )/16)j , 

V 2fi = - Sm( ^ 2) exp(-100(x 2 + x 2 )) exp (l(- Xl + cos(4a;2)/16)) . 

The final time is T — 1.0. We take e — 1/64. Figure [3] compares the pressure II 
of the true solution with the one by EFGA. Figure |4] compares the velocity V of 
the true solution with the one by EFGA. It is clear that EFGA can provide a good 
approximation to both the pressure and velocity for acoustic wave propagation in 
two dimension. The true solution is given by the spectral method using the mesh 
Sxi = 6x2 = 1/512 for domain [—1.5,0.5] x [—1,1]. We take Sqi = 8q 2 = 8p\ = 
5p 2 = S yi = 6y 2 = 1/32 and 5x x = 5x 2 = 1/64 in EFGA. 

5.2. Schrodinger equation. 

Example 5.3 (One-dimensional Schrodinger equation). 

e 2 



2(5 
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Figure 2. Example the comparison of the true solution (solid 
line) and the solution by EFGA (dash line). Left: the real part of 
wave field; right: the errors between them. 



and the initial condition is 



exp 



x 
2s 



We use this one-dimensional Schrodinger equation with zero potential as an 
example to compare the performance of Lagrangian and Eulerian methods. The 
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Figure 3. Example I5.2[ the comparison of the true solution and 
the solution by EFGA. Top (left): pressure of EFGA; top (right): 
pressure of true solution; bottom: the error between them. 

true solution can be given analytically, 

*'<M>^^H + 4«)> 

which implies the solution spreads as time increases. 

We choose e = 1/128 and evolve the equation up to T = 10. The mesh sizes 
are 5x = Sy = Sq = Sp = 1/64. We take N q — 64 and N p = 33 in the Lagrangian 
method. The comparison of wave amplitudes and numerical errors are presented 
in Figure [5j One can see that, when the divergence of particle trajectories occurs, 
Eulerian method has a much better resolution than the Lagrangian method. 

Example 5.4 (Two-dimensional Schrodinger equation). 

and the initial condition is 

= exp(— 25(x\ + x%)) exp sin(xi) sin^)) . 

This is an example of Schrodinger equation in two dimension, which describes 
the dynamics of electron under harmonic potential. We take e = 1/128. Figure [S] 
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1 1 

(a) Eulerian Frozen Gaussian approximation 




1 -0.8 -0.6 -0.4 -0.2 



(c) Errors 

Figure 4. Example I5.2[ the comparison of the true solution and 
the solution by EFGA. Left: velocity component V\] right: velocity 
component V%. 



compares the wave amplitude of the true solution with the numerical one at time 
T = 0.5 and T = 1. This shows that Eulerian Herman-Kluk propagator has good 
performances in both cases of solution spreading and localizing. The true solution 
is given by the spectral method using the mesh Sx\ = 8x2 = 1/512 for domain 
[—1,1] x [—1,1]. In the numerical approximation, the mesh sizes are chosen to be 
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Figure 5. Example 15.31 Left: the wave amplitude comparison of 
the true solution (dashed line), the Lagrangian method (dots) and 
the Eulerian method (solid line); right: the numerical errors of the 
Lagrangian method (dots) and the Eulerian method (solid line). 

Sqi = Sq2 = 5pi = 5p2 = 5y± = 5y2 = 1/32 in discretization of integrals and 
5xi = 5x2 = 1/32 in reconstruction of solution. 

6. Conclusion 

We extend the formulation of frozen Gaussian approximation to general linear 
strictly hyperbolic system. Based on the Eulerian formulation of frozen Gaussian 
approximation, Eulerian methods are developed to resolve the divergence problem 
of the Lagrangian method. Moreover, the Eulerian methods can be also used for 
computing the Herman-Kluk propagator of the Schrodinger equation in quantum 
mechanics. The performance of the proposed methods is verified by numerical 
examples. This paper, together with [15] . provides an efficient methodology for 
computing high frequency wave propagation for general hyperbolic systems with 
smooth coefficient. 
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(a) Eulerian Herman-Kluk propagator 
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(b) True solution 
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Figure 6. Example I5.4[ the comparison of the true solution and 
the solution by Eulerian Herman-Kluk propagator. Left: T = 0.5; 
right: T = 1. 
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